R包vegan的冗余分析（RDA）
某研究对三种环境中的土壤进行了采样（环境A、B、C，各环境下采集9个土壤样本），结合二代测序观测了土壤细菌群落物种组成，并对多种土壤理化指标进行了测量。获得两个数据集。
（1）文件“phylum_table.txt”为通过16S测序所得的细菌类群丰度表格，此处统计至细菌“门水平”（即界门纲目科属种的“门”这一水平）。（为了方便演示，我没使用OTU水平的）
（2）文件“env_table.txt”记录了多种土壤理化指标信息，即环境数据。
接下来期望通过RDA分析这些数据，用以分析群落物种组成与环境梯度的关系。
RDA的初步执行
vegan包是生态学统计分析中最常用的R包之一，已经打包好了关于RDA计算、检验等一些列命令，方便我们快速有效地执行RDA。
vegan包中，通过rda()执行RDA，该函数的输入格式大致如下。
library(vegan)#详情 ?rda
#调用格式 1
rda(Y, X, W)
#或者格式 2
rda(Y~var1+var2+var3+factorA+var2*var3+Condition(var4))
在第1种调用格式中，Y为响应变量矩阵，X为解释变量矩阵，W为偏RDA分析需要输入的协变量矩阵。当不考虑协变量时，可直接输入响应变量矩阵Y和解释变量矩阵X，即“rda(Y, X)”。这种写法虽然简单，但要求解释变量矩阵或协变量矩阵中不能含有因子变量（定性变量）。（当仅存在响应变量矩阵Y时，则执行PCA排序）
因此，通常建议使用第2种写法。其中Y为响应变量矩阵，var1、var2、var3为定量解释变量，factorA为因子变量，var2*var3为变量2和变量3的交互作用项，var4为协变量。与上述第1种写法相比，不仅能够识别因子变量，还能够对变量间的交互作用加以考虑。下文的具体示例中，均以第2种写法作展示。
RDA包含RDA、基于转化的冗余分析（tb-RDA）、偏冗余分析（偏RDA）、基于距离的冗余分析（db-RDA）、非线性RDA等多种模式。实际的数据分析中，需根据情况选择合适的分析模式。关于这些RDA的基本概念描述，可参考前文。
以下分别简介在R中运行这几种RDA模式的命令方法。
常规RDA
即提供响应变量矩阵Y，以及解释变量矩阵X，将原始数据直接通过rda()函数分析。
##读取数据
#读取物种数据，细菌门水平丰度表
phylum <- read.delim('phylum_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#读取环境变量
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#使用全部的环境数据
rda_result <- rda(phylum~., env, scale = FALSE)
使用物种数据和环境变量数据，其中原始物种组成数据（数据框phylum）作为响应变量，原始环境因子数据（数据框env）作为解释变量。式中scale参数用于对响应变量矩阵（与解释变量矩阵无关）执行标准化处理。scale = FALSE，不执行标准化，使用拟合值的协方差矩阵运算RDA；若为TRUE，则执行标准化，基于相关矩阵运算RDA。响应变量的量纲不一致时（不过，物种多度数据不存在这种情况），或者存在极端值影响时（也可以通过额外的数据转化来解决），推荐scale = TRUE。
“~.”表示将数据框env中所有的列变量（环境因子）作为解释变量带入RDA排序，是一种变量全选的省略模式（在不考虑协变量以及变量间的交互作用项的情况下，可以直接使用），能够自动识别定量解释变量以及定性解释变量类型。对于特殊情况，可以写作如下样式。
#若只关注局部环境数据，除了在原始表格中修改变量个数外，还可直接在 rda() 中指定#例如只考虑 pH、TC、TN、AP、AK 这 5 种环境变量
rda_result <- rda(phylum~pH+TC+TN+AP+AK, env, scale = FALSE)
#默认情况下，rda(phylum~., env)，不包含环境变量间的交互作用
#若想关注某两个环境变量间的交互作用，例如添加 TC 和 TN 的交互作用项
rda_result <- rda(phylum~pH+TC+DOC+SOM+TN+NO3+NH4+AP+AK+TC*TN, env, scale = FALSE)
#当存在协变量时（即偏 RDA，可参见下文）
#例如控制土壤 pH 影响后（pH 作为协变量），观测其它环境因素的影响
rda_result <- rda(phylum~TC+DOC+SOM+TN+NO3+NH4+AP+AK+Condition(pH), data = env, scale = FALSE)
此时我们将RDA结果存储在对象rda_result中，内容说明参见下文“RDA结果的初步解读”。
tb-RDA（本篇重点展示的是它）
对应于tb-PCA，通常情况下，对于物种多度数据，不建议直接用于RDA排序。特别是当物种丰度表中出现很多“0”值的情况下，应当首先进行Hellinger转化后，再输入至RDA中。这种首先对原始数据作一定的转化后在执行RDA的方法称为基于转化的冗余分析（Transformation-based redundancy analysis，tb-RDA）。除了第一步的数据转化外，其与过程均和常规的RDA相同。
因此我们首先对示例物种丰度数据进行Hellinger转化，并使用转化后的数据执行RDA。
##tb-RDA
#物种数据 Hellinger 预转化（处理包含很多 0 值的群落物种数据时，推荐使用）
phylum_hel <- decostand(phylum, method = 'hellinger')
#使用全部的环境数据
rda_tb <- rda(phylum_hel~., env, scale = FALSE)
#参考下文中对 RDA 结果解读的说明，可分别使用 plot(rda_result) 和 plot(rda_tb) 查看比较 Hellinger 转化前后二者的差异
rda()命令中，使用Hellinger转化后的物种数据（数据框phylum_hel）作为响应变量，原始环境因子数据（数据框env）作为解释变量。式中scale参数用于对响应变量矩阵（与解释变量矩阵无关）执行标准化处理。scale = FALSE，不执行标准化，使用拟合值的协方差矩阵运算RDA；若为TRUE，则执行标准化，基于相关矩阵运算RDA。Hellinger转化后的物种多度数据，通常无需标准化，首先物种数据的量纲均相同，其次Hellinger转化消除了极端值的影响，降低了高丰度物种的权重。
此时我们将tb-RDA结果存储在对象rda_tb中，内容说明参见下文“RDA结果的初步解读”。
偏RDA
有时，我们期望在控制变量矩阵W对响应变量Y的影响的前提下，将关注点集中在另一组变量矩阵X对响应变量Y的影响上，应用到RDA中即为偏冗余分析（Partial canonical ordination ，偏RDA）。
例如对于示例数据，我们期望控制土壤pH影响后（pH作为协变量），观测其它环境因素的影响。
##偏 RDA
#控制土壤 pH 影响后（pH 作为协变量），观测其它环境因素的影响；物种数据 Hellinger 预转化
rda_part <- rda(phylum_hel~TC+DOC+SOM+TN+NO3+NH4+AP+AK+Condition(pH), data = env, scale = FALSE)
rda()命令中，使用Hellinger转化后的物种数据（数据框phylum_hel）作为响应变量，原始环境因子数据（数据框env）作为解释变量。由于加入了协变量，此时“~”后不再直接使用“.”表示全部的环境因子，而是通过指定的方式，将所考虑的解释变量的名称以“+”连接，并将协变量（示例为pH）添加至Condition()里。
此时我们将偏RDA结果存储在对象rda_part中，内容说明参见下文“RDA结果的初步解读”。
在变量选择过程中，经常使用到偏RDA，变量选择的具体示例详见下文。
db-RDA
基于距离的冗余分析（Distance-based redundancy analysis，db-RDA）。
考虑到它比较特殊，且又非常的重要，本篇暂且跳过这种模式，我把它单独整理了一篇，放在下篇展示。
非线性关系的RDA
当响应变量与解释变量明显为非线性关系（例如单峰关系，当然，若响应变量与解释变量存在明显的单峰响应模式，可使用CCA来解决），常规的线性RDA模型不适合时，可考虑使用这种方法来解决。
由于我并没有充分理由认为示例数据中响应变量与解释变量存在明显的非线性关系（鉴定过程繁琐，懒得观测……），故这里不再使用示例数据演示这种非线性关系的RDA过程了。大家有兴趣可参考赖江山老师译作“数量生态学:R语言的应用”170-174页内容。
RDA结果的初步解读和信息提取
上文简介了使用R的vegan包中的rda()命令，初步获得几种不同模式RDA排序结果的方法。在初步得到排序结果后，我们首先需要对结果进行解读。
下文仅对rda()运行结果中的各部分内容所反映的信息（意义）及其提取方法作简要描述，并初步展示排序图的绘制方法，即主要以软件所得运行结果为例展开叙述以帮助快速上手软件使用。对于RDA排序结果解读的更详细说明，即主要的基本概念部分，如变差（方差）解释程度、模型初步评估、排序图的详解、I型或II型标尺的缩放原理等，本文不作介绍，可参考前文。
RDA结果总概
以上述tb-RDA排序结果“rda_tb”为例，对vegan所得RDA结果的主要结构内容进行简要的说明。
##RDA 结果解读，以下以 tb-RDA 结果为例#查看统计结果信息，以 I 型标尺为例
rda_tb.scaling1 <- summary(rda_tb, scaling = 1)
rda_tb.scaling1
summary()展示结果时，根据所重要要关注的内容，可选根据I型标尺或II型标尺缩放排序坐标后展示出。如果对排序样方之间的距离更感兴趣，或者大多数解释变量为因子类型，则考虑I型标尺；如果对变量之间的相关关系更感兴趣，则考虑II型标尺。关于I型标尺和II型标尺，可参考前文。这里我们展示了按I型标尺缩放后排序结果（参数scaling = 1；同理，scaling = 2展示II型标尺；参数中还有scaling = 0和scaling = 3的选项，这个我就不清楚了），部分内容如下所示。
仅通过初步结果，我们发现该RDA模型承载的总方差比例达到了70%以上（0.7338），特别是RDA1和RDA2的解释量总和为“0.53445+0.118470=0.65292”（即承载了响应变量矩阵65.29%的方差），非常可观的一个数字。做到这里是不是很开心地去整理数据去了？别急，仔细看下红字部分的标注，这些结果并不是很可信，有待于校正和检验通过后才能被接受。详见下文“RDA的R2校正和约束轴的置换检验”。
对于偏RDA结果，解读方式也大致同上。但是有这么几个地方存在区别。
RDA排序图
先不急介绍RDA模型的验证方法，还得先把基础部分介绍完…...
排序结果最终以排序图作为呈现。对于vegan所得rda排序结果，可直接使用plot()绘制排序图。若想提前了解更多，可使用?plot.cca()查看详细说明（注意不是?plot()）。
此处同样以上述tb-RDA排序结果“rda_tb”为例，绘制简单样式的RDA三序图作为展示。图中对于排序对象（样方），展示为“标签点”；响应变量（物种变量）和解释变量（环境变量）均以向量表示。
#作图查看排序结果，三序图，包含 I 型标尺和 II 型标尺par(mfrow = c(1, 2))
plot(rda_tb, scaling = 1, main = 'I 型标尺', display = c('wa', 'sp', 'cn'))
rda_sp.scaling1 <- scores(rda_tb, choices = 1:2, scaling = 1, display = 'sp')
arrows(0, 0, rda_sp.scaling1[ ,1], rda_sp.scaling1[ ,2], length = 0, lty = 1, col = 'red')
plot(rda_tb, scaling = 2, main = 'II 型标尺', display = c('wa', 'sp', 'cn'))
rda_sp.scaling2 <- scores(rda_tb, choices = 1:2, scaling = 2, display = 'sp')
arrows(0, 0, rda_sp.scaling2[ ,1], rda_sp.scaling2[ ,2], length = 0, lty = 1, col = 'red')
plot()命令中，display参数用于定义在排序图中展示哪些信息。其中，“sp”代表物种，“wa”代表使用物种加权和计算的样方坐标，“lc”代表以拟合值（解释变量的线性组合）计算样方的坐标，“cn”代表约束成分（即解释变量）。可知对于样方坐标有两种可选展示方式，上述示例使用了“wa”。选择哪种排序坐标在排序图绘制样方的点，至今还有争议（参见“数量生态学:R语言的应用”148页内容），但无论如何，排序图的解读都需要依据统计显著性检验结果，与多元回归一样，不显著的回归结果不能被解读，必须丢弃（RDA的统计检验部分参见下文）。
plot()默认作图结果中，物种同样以“标签点”来表示，因此我们需要额外绘制箭头以展示物种向量。scores()用于提取物种、样方或环境变量的排序坐标，同样以display参数定义。这里使用scores()提取了物种变量的排序坐标，并使用arrows()在排序图中为物种向量添加箭头，以(0,0)出发，指向对应的坐标位置处。
我们简单地绘制了I型标尺和II型标尺的RDA三序图，结果如下所示。关于排序图的解读，可参考前文。
由于物种信息较多，影响观测。如果仅考虑通过排序图展示环境变量对样方的约束信息，那么通常可以在图中忽略物种信息，即展示为双序图样式。
下图展示了仅包含样方和环境变量信息的双序图，同时分别使用display = c('wa', 'cn')和display = c('lc', 'cn')参数，在排序图中分别展示了使用物种加权和计算的样方坐标以及以拟合值（解释变量的线性组合）计算样方的坐标。（可观察比较二者的区别）
#隐藏物种信息，以 I 型标尺为例展示双序图，并查看分别使用物种加权计算的样方坐标（wa）以及拟合的样方坐标（lc）的差异par(mfrow = c(1, 2))
plot(rda_tb, scaling = 1, main = 'I 型标尺，加权', display = c('wa', 'cn'))
plot(rda_tb, scaling = 1, main = 'I 型标尺，拟合', display = c('lc', 'cn'))
以上暂时仅对排序图的可视化做个简单介绍。对于更详细的可视化细节部分，可参见本篇末尾。
备注：对于vegan的排序结果使用plot()绘制排序图时，所得排序图中解释变量（环境变量）向量的箭头顶点坐标可能并非我们通过summary()所观测到的实际坐标位置，而更像是在原坐标的基础上同比例放大（或缩小）了一定的数值后展示的。这是因为在某些情况下，若直接展示原始坐标，则向量箭头的长度相较于排序图将会显得比例失调（过长或过短），严重影响我们对图的解读，此时通常考虑放大（或缩小）一定数值后展示出（可参见Legendre和Legendre，1998，404页）。vegan默认的做图方式中即考虑了这一点。
RDA结果中的主要信息提取
若有需要，考虑在结果中将需要的信息提取出来。以下同样以上述tb-RDA排序结果“rda_tb”为例展示RDA重要信息的提取方法。
上述在初步绘制排序图观测时，提到了可使用scores()提取物种、样方或环境变量的排序坐标，并且已经使用scores()提取了物种变量的排序坐标。同样地，修改display参数后即可提取样方和环境变量的排序坐标。
#scores() 提取排序得分（坐标），以 I 型标尺为例，前四轴为例#使用物种加权和计算的样方得分
rda_site.scaling1 <- scores(rda_tb, choices = 1:4, scaling = 1, display = 'wa')
#物种变量（响应变量）得分
rda_sp.scaling1 <- scores(rda_tb, choices = 1:4, scaling = 1, display = 'sp')
#环境变量（解释变量）得分
rda_env.scaling1 <- scores(rda_tb, choices = 1:4, scaling = 1, display = 'bp')
#或者在 summary() 后提取，以 I 型标尺为例，前四轴为例
rda_tb.scaling1 <- summary(rda_tb, scaling = 1)
names(rda_tb.scaling1)
#使用物种加权和计算的样方得分
rda_site.scaling1 <- rda_tb.scaling1$site[ ,1:4]
#物种
rda_sp.scaling1 <- rda_tb.scaling1$species[ ,1:4]
#环境
rda_env.scaling1 <- rda_tb.scaling1$biplot[ ,1:4]
#若需要输出在本地
#样方
write.table(data.frame(rda_site.scaling1), 'rda_site.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#物种
write.table(data.frame(rda_sp.scaling1), 'rda_sp.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#环境
write.table(data.frame(rda_env.scaling1), 'rda_env.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#不建议直接在原始数据集中提取，因为这样提取的坐标数值未经标尺缩放处理，不利于反映生物学问题
#rda_tb$CCA$u[ ,1:4]
#rda_tb$CCA$v[ ,1:4]
#rda_tb$CCA$biplot[ ,1:4]
RDA基于多元回归，若对典范特征系数（即每个解释变量与每个约束轴之间的回归系数）感兴趣，可使用coef()提取（使用?coef()以了解更多）。
#coef() 提取典范系数rda_coef <- coef(rda_tb)
关于约束模型解释的总变差（即响应变量矩阵的总方差能够被解释变量解释的部分，如果用比例表示，即等同于多元回归的R2，也就是约束轴的总解释量）以及各约束轴解释量等，不建议直接在原始结果中获取。原因很简单，正如前文所提到，初始所得的RDA模型需要经过校正和检验后才可使用，否则可能会带来较大的偏差。
因此，下文的内容很重要，对于严谨的RDA模型的构建，请不要忽略。
RDA的R2校正和约束轴的置换检验
下文仅对vegan实现RDA的R2校正以及约束轴的置换检验的方法（函数命令）作简介，基本概念性质的内容，可参考前文。
R2校正
上文中提到，与多元回归的未校正R2一样，RDA的R2也有偏差，需要进行校正（即原始R2不可直接使用）。RDA结果中所展示的各约束轴的解释量也是基于未校正的R2得来的，也需要根据校正后的R2重新计算（即原始所得约束轴的解释量是不可信的，同样不能直接使用）。
对于校正前后R2的提取，可使用RsquareAdj()完成。
#RsquareAdj() 提取 R2，详情 ?RsquareAdj()
r2 <- RsquareAdj(rda_tb)
rda_noadj <- r2$r.squared #原始 R2
rda_adj <- r2$adj.r.squared #校正后的 R2
R2经过校正后通常会降低，这是必然的。
对于示例数据的tb-RDA排序结果，通过上文，我们已经得知了响应变量矩阵的总方差为0.02686，未校正前的约束轴解释的方差总和为0.01971（占比0.7338，即未校正前的R2）；那么根据校正后的R2（0.5928），可还原校正R2后的约束轴解释的方差总和约为“0.02686×0.5928=0.01592”。
同样地，我们也可获得各约束轴校正后的解释量。例如，已知RDA第一个约束轴（按解释量或特征值的大小由高往低排序，各约束轴记为RDA1、RDA2、RDA3……）的解释方差占约束模型总解释方差的0.72835（特征值占比0.72835，同样地可知，该轴的解释量占约束模型总解释量的比例也是0.72835），由于校正R2后的约束模型的总解释量降至0.5928，那么可推断校正R2后的第一个约束轴的解释量（第一轴的解释方差占响应变量矩阵总方差的比例）即为“0.5928×0.72835=0.43177”。同理，校正R2后RDA第二个约束轴的解释量“0.5928×0.16145=0.09571”，其余约束轴的解释量以此类推。
对于我们的示例数据的tb-RDA结果，尽管在R2校正后约束模型的总体解释量有所下降（由0.7338下降至0.5928），但看起来似乎仍很可观。特别是的前两个约束轴的解释量仍然非常好：RDA1和RDA2的解释量总和为“0.43177+0.09571=0.52748”，即承载了响应变量矩阵52.75%的方差，达到了一半以上。
好了，这时候，我们的RDA模型可靠了吗？不，还需检验约束轴的显著性。
约束轴的置换检验及p值校正
首先需对全模型执行置换检验查看显著性，仅当全模型通过检验时，才能依次对每一个约束轴执行检验判断是否能够保留。否则，RDA模型则从全局上来讲就不合理，本次的RDA结果不能再被使用，需要考虑更改分析方案（更换排序模型，或者再尝试对数据进行有效的转化后使用等）。
vegan包中执行置换检验的命令是anova()（很尴尬它就叫这个名字，R语言中方差分析的命令也是anova()……二者名称一样但用法完全不同，注意区分二者的用途）。关于此命令的详情，可使用?anova.cca()查看帮助（不要使用?anova()，不然弹出的将会是方差分析的命令帮助）。
以下以前文tb-RDA结果为例，展示RDA全模型及各约束轴的置换检验过程。
##置换检验
#所有约束轴的置换检验，基于 999 次置换检验
rda_tb_test <- anova(rda_tb, permutations = 999)
#或者使用
rda_tb_test <- anova.cca(rda_tb, step = 1000)
#各约束轴逐一检验，基于 999 次置换检验
rda_tb_test_axis <- anova(rda_tb, by = 'axis', permutations = 999)
#或者使用
rda_tb_test_axis <- anova.cca(rda_tb, by = 'axis', step = 1000)
若使用anova()，使用permutations参数指定置换次数；若使用anova.cca()，使用step参数设定产生参照分布的样方总数，step数值减1即为实际的置换次数。对于RDA全模型的置换检验，直接将前文所得tb-RDA结果“rda_tb”输入至命令内即可；对于各约束轴，则还需在命令中添加by = 'axis'。此外，anova()还可以通过by = ' margin'检验协变量等，这里不再多说，有需要请?anova.cca()参阅该命令的帮助文档。
结果就很明了了。首先全模型是显著的，因此我们继续检验每个约束轴，发现仅第一轴（RDA1）和第二轴（RDA2）是显著的（默认显著性阈值p<0.05），因此第三轴及以后的结果需要被剔除。
当然，必要时可以考虑做下p值校正。
对于上述示例结果，我们加入Bonferroni校正，使用p.adjust()完成。
#p 值校正（Bonferroni 为例）rda_tb_test_axis$`Pr(>F)` <- p.adjust(rda_tb_test_axis$`Pr(>F)`, method = 'bonferroni')
结果如下所示。对于RDA1和RDA2两个约束轴，结果仍然是可信的。
对于示例数据来讲，根据检验结果，我们仅能在全模型中提取前两个约束轴的信息出来以进行更进一步的分析。不过，这样的结果似乎已经是挺不错的了。如上文所述，在校正了约束模型的R2后，RDA1和RDA2的解释量总和为“0.43177+0.09571=0.52748”（即承载了响应变量矩阵52.75%的方差，达到了一半以上；再结合置换检验结果，即便不考虑其它的约束轴，约束模型所能解释的方差占比也是蛮可观的），表明给定的环境变量对群落物种组成分布具有相当高的解释度。特别是第一轴，提示我们在解读排序图时，可能需要主要观察样方或物种沿第一约束轴的分布规律，以及相关环境变量与约束轴、样方及物种的相关性程度等。
事实上，即便后面的约束轴通过了检验，我们一般也不将它们列为考虑的范畴，毕竟它们的解释量很低，强制将它们加入生态学模型解释物种组成成因也可能会带来混乱的结果。我们更期望的还是具有较高解释量的约束轴能够被保留，以及样方或物种的排序分布具有明显趋势的模型不要被舍弃等。
注意，非约束模型中排序轴取舍的常用方法不适用于约束模型。例如断棍模型、凯撒格特曼规则（Kaiser-Guttman rule）等，我们可能在PCA分析中常用它们作为排序轴的选择标准，但它们在RDA约束轴的取舍中意义不大。对于约束轴的取舍，还需通过置换检验获得。然而，对于RDA的非约束残差轴，断棍模型等仍然适用。
非约束残差轴的重要性确定
当我们发现RDA模型的承载方差比例较低，存在相当一部分比例的残差未参与解释时，可能需要审视非约束残差轴。以下接着以上文tb-RDA结果为例，提取其中非约束残差轴的特征值，并分别依据断棍模型和Kaiser-Guttman准则确定残差轴。
##断棍模型和 Kaiser-Guttman 准则确定残差轴
#提取残差特征值
pca_eig <- rda_tb$CA$eig
#Kaiser-Guttman 准则
pca_eig[pca_eig > mean(pca_eig)]
#断棍模型
n <- length(pca_eig)
bsm <- data.frame(j=seq(1:n), p = 0)
bsm$p[1] <- 1/n
for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i))
bsm$p <- 100*bsm$p/n
bsm
# 绘制每轴的特征根和方差百分比
par(mfrow = c(2, 1))
barplot(pca_eig, main = '特征根', col = 'bisque', las = 2)
abline(h = mean(pca_eig), col = 'red')
legend('topright', '平均特征根', lwd = 1, col = 2, bty = 'n')
barplot(t(cbind(100 * pca_eig/sum(pca_eig), bsm$p[n:1])), beside = TRUE, main = '% 变差', col = c('bisque', 2), las = 2)
legend('topright', c('% 特征根', '断棍模型'), pch = 15, col = c('bisque', 2), bty = 'n')
结果如下。综合考虑，大约3个残差非约束轴是比较重要的。若有必要，我们可能需要绘制前2-3个非约束轴的排序图，查看为何还有相当一部分方差未能解释。
变量选择（前向选择）
约束排序中，有时解释变量太多，需要想办法减少解释变量的数量。一般来讲，减少解释变量的数量可能有两个并不冲突的原因：第一是寻求简约的模型，利于我们对模型的解读；第二是当解释变量过多时会导致模型混乱，例如有些解释变量之间可能存在较强的线性相关，即共线性问题，可能会造成回归系数不稳定。
下文仅简要展示变量选择过程在R中的实现方法，关于RDA变量选择的更详细原理和概念描述等，可参考前文。
方差膨胀因子
每个变量的共线性程度可以使用变量的方差膨胀因子（Variance Inflation Factor，VIF）度量，VIF是衡量一个解释变量的回归系数的方差由共线性引起的膨胀比例。如果VIF值超过20，表示共线性很严重。实际上，VIF值超过10则可能就会有共线性的问题，需要处理。
vegan包中计算约束模型中变量VIF值的命令为vif.cca()。以上文tb-RDA结果为例计算VIF。
#计算方差膨胀因子
vif.cca(rda_tb)
可以看出有几个变量的VIF值特别大（超过10甚至20），所以有必要剔除一些变量。
多元回归变量筛选通常有三种模式：前向选择（forward selection）、后向选择（backward selection）以及逐步选择（stepwise selection，也称双向选择，forward-backward selection）。其中前向选择在RDA分析中最为常用，以下将以前向选择为例简介RDA模型中的变量选择方法。
前向选择
在R中，解释变量的前向选择通常使用vegan包中ordistep()函数、ordiR2step()函数，或者packfor包中的forward.sel()函数来完成。不同的命令执行前向选择所依据的原理是不同的，以下分别简介用法。
ordistep()
vegan包中ordistep()函数执行前向选择过程大致如下（除了前向选择，ordistep()同样可以执行后向选择和逐步选择，本篇不再细说，详细请使用?ordistep()查看说明）。
（1）依次分别运行每个解释变量与响应变量的RDA排序。
（2）选择“最好”的解释变量。对于ordistep()，根据置换检验中F统计量的显著性水平（p值）选择最好的变量，即最小的p值。如果遇到p值相等的情况，拥有最小AIC信息准则（Akaike information criterion，AIC）的变量应该入选。最好的变量也就是最显著的变量。
（3）寻找模型中第二个、第三个、第四个……解释变量。上一步选取了只含有一个最好变量的模型，现在再加入某一个剩下的解释变量。对于ordistep()，选择偏RDA置换检验最小p值的模型，即表示新加入的变量对模型的贡献最显著。如果遇到p值相等的情况，具有最小AIC值的模型应该入选。
（4）这个过程一直持续，直到无显著的解释变量为止。即如果加入新变量的偏RDA置换检验显著性p值大于或等于α（例如，以α=0.05为准），选择过程即被终止。
在上文，通过计算计算方差膨胀因子，我们已知tb-RDA示例结果中的某些解释变量最好加以剔除。以下使用ordistep()，展示在示例tb-RDA模型中执行前向选择的过程。ordistep()可以直接将rda()的输出作为输入（此外，ordistep()也同样可用于CCA模型，vegan包cca()的输出）。由于原始的tb-RDA模型就已经是包含所有解释变量的RDA全模型了，我们可直接在它基础上继续执行前向选择。
#vegan 包 ordistep() 前向选择，基于 999 次置换检验，详情 ?ordistep()rda_tb_forward_p <- ordistep(rda(phylum_hel~1, env, scale = FALSE), scope = formula(rda_tb), direction = 'forward', permutations = 999)
#简要绘制双序图比较变量选择前后结果
par(mfrow = c(1, 2))
plot(rda_tb, scaling = 1, main = '原始模型，I 型标尺', display = c('wa', 'cn'))
plot(rda_tb_forward_p, scaling = 1, main = '前向选择后，I 型标尺', display = c('wa', 'cn'))
#细节部分查看
summary(rda_tb_forward_p, scaling = 1)
#比较选择前后校正后 R2 的差异
RsquareAdj(rda_tb)$adj.r.squared
RsquareAdj(rda_tb_forward_p)$adj.r.squared
ordistep()中，rda(phylum_hel~1, env)意为从env（环境解释变量数据集矩阵）中的第一个变量开始执行选择；direction = 'forward'，前向选择；permutations = 999，999次置换。
全模型中共计9个环境变量，经过ordistep()前向选择后保留了其中的6个。尽管解释变量的数量少了1/3，但是我们发现模型的总解释量几乎没有改变，即变量选择后在提供了简约的模型的同时，并没有牺牲解释能力。
然后再根据排序图，或者summary()统计，重新解读模型，观察比较变量选择前后的差异；并重新校正模型R2以及检验约束轴，必要时增添p值校正，最终选择合适的约束轴用于解释生态学现象等。此处不再更多的展示，请大家自行测试。
但是ordistep()的选择标准过于宽松，有时会选择显著但不包含任何变量的模型（夸大I类错误），或选择包括过多变量的模型（夸大被解释方差的量）。为了解决这个问题，可将Borcard终止准则（Blanchet等，2008）引入作为前向选择的第二个终止原则，加强版的前向选择程序，也就是下文将要介绍的ordiR2step()和forward.sel()。
在使用ordistep()时，若是发现结果似乎异常，例如下图中的这种结果，变量选择后的校正R2高于了全模型的校正R2？则推荐使用ordiR2step()或forward.sel()执行变量选择过程（其实也推荐在一开始就直接使用这种加强版的前向选择程序ordiR2step()或forward.sel()，忽略ordistep()）。
ordiR2step()
ordiR2step()的前向选择原理与ordistep()类似，但引入全模型R2adj作为第二个终止原则（即增添了Borcard终止准则）：如果当前所选模型的R2adj达到或超过全模型的R2adj，或备选变量的偏RDA置换检验p值不显著，则变量选择将停止。
与上文一致，以下仍然以示例tb-RDA模型为例做演示。和ordistep()用法一致，除了前向选择，ordiR2step()也可以执行后向选择和逐步选择，并且可以直接将rda()的输出作为输入（此外，ordiR2step()也同样可用于CCA模型，vegan包cca()的输出）。本文不再细说，详细请使用?ordiR2step()查看说明。
#vegan 包 ordiR2step() 前向选择，基于 999 次置换检验，详情 ?ordiR2step()rda_tb_forward_r <- ordiR2step(rda(phylum_hel~1, env, scale = FALSE), scope = formula(rda_tb), R2scope = rda_adj, direction = 'forward', permutations = 999)
#简要绘制双序图比较变量选择前后结果
par(mfrow = c(1, 2))
plot(rda_tb, scaling = 1, main = '原始模型，I 型标尺', display = c('wa', 'cn'))
plot(rda_tb_forward_r, scaling = 1, main = '前向选择后，I 型标尺', display = c('wa', 'cn'))
#细节部分查看
summary(rda_tb_forward_r, scaling = 1)
#比较选择前后校正后 R2 的差异
RsquareAdj(rda_tb)$adj.r.squared
RsquareAdj(rda_tb_forward_r)$adj.r.squared
ordiR2step()中，rda(phylum_hel~1, env)意为从env（环境解释变量数据集矩阵）中的第一个变量开始执行选择；rda_adj即为RDA全模型的校正R2，在上文中已经通过RsquareAdj()获得，参数scope = formula(rda_tb)即指定全模型的校正R2作为判断的依据；direction = 'forward'，前向选择；permutations = 999，999次置换。
全模型中共计9个环境变量，经过ordistep()前向选择后保留了其中的6个。尽管解释变量的数量少了1/3，但是我们发现模型的总解释量几乎没有改变，即变量选择后在提供了简约的模型的同时，并没有牺牲解释能力。
然后再根据排序图，或者summary统计，重新解读模型，观察比较变量选择前后的差异；并重新校正模型R2以及检验约束轴，必要时增添p值校正，最终选择合适的约束轴用于解释生态学现象等。此处不再更多的展示，请大家自行测试。
备注：在本示例中，ordiR2step()和ordistep()的所得结果是相同的，巧合了而已（本身数据集不大，解释变量的数量也不多）。但在实际的数据分析中，二者的结果有时差别还是挺大的。
提取变量选择后的排序坐标，以此为例。
#提取或输出变量选择后的排序坐标#以 ordiR2step() 结果为例，以 I 型标尺为例，前四轴为例
#提取方式可参考上文，这里通过 scores() 提取
#使用物种加权和计算的样方得分
rda_tb_forward_r_site.scaling1 <- scores(rda_tb_forward_r, choices = 1:4, scaling = 1, display = 'wa')
write.table(data.frame(rda_tb_forward_r_site.scaling1), 'rda_tb_forward_r_site.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#物种变量（响应变量）得分
rda_tb_forward_r_sp.scaling1 <- scores(rda_tb_forward_r, choices = 1:4, scaling = 1, display = 'sp')
write.table(data.frame(rda_tb_forward_r_sp.scaling1), 'rda_tb_forward_r_sp.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#环境变量（解释变量）得分
rda_tb_forward_r_env.scaling1 <- scores(rda_tb_forward_r, choices = 1:4, scaling = 1, display = 'bp')
write.table(data.frame(rda_tb_forward_r_env.scaling1), 'rda_tb_forward_r_env.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
forward.sel()
由于forward.sel()和ordiR2step()的选择原理相似，都是引入了Borcard终止准则后的加强版的选择程序，所以一般情况下二者前向选择的结果是一致的（尽管具有不同的参数设置逻辑）。
但是与ordiR2step()相比，forward.sel()目前仅适用于RDA（改良后的ordiR2step()既可用于RDA，vegan包rda()的输出；也可用于CCA，vegan包cca()的输出），且仅可用前向选择（ordiR2step()同时支持前向、后向、逐步选择），所以看起来似乎并没有ordiR2step()好用……
而且forward.sel()所在的packfor包仅支持特定的R版本，所以本篇就不展示它了，大家有兴趣可自行研究。
变量选择中的注意事项
尽管变量选择策略具有明显的优势，但一定切记不应盲目依赖自动选择程序在回归模型中选择相关的环境变量，因为可能得到生态学上不相关的模型，或者被忽视的其它变量组合（当前保留变量集中的部分解释变量+某些被剔除的解释变量）在某种程度上也可能比当前模型更具代表性（Legendre和Legendre，1998）。即仅通过统计学手段得到的最优变量集合可能并非代表了最具生态学意义的模型。
排序图的可视化示例
到这里，对RDA分析的介绍基本就结束了。其实到了这一步，我们期间通过R2校正、置换检验、前向选择步骤等，再结合自己的专业知识进行合理的解读和评估后，也差不多得到了一个合理的RDA模型了。当然，我的意思可不是指这里的示例数据的RDA分析已经得到了一个理想模型了，而是指在实际情况的分析中，我们通过多方考虑所得的最终模型，这个模型除了在统计学上合理外，在解释生物学意义中也是较为理想的。这个最终的模型和最初直接由rda()所得模型相比，可能差异巨大。不过到了这里，想必大家也能理解了为什么最初所得的RDA模型最好不要直接使用的原因了。
那么最后，再对RDA排序图的可视化方案做个简介。
plot()作图
在前文，已经简单地使用plot()函数（实为plot.cca()）绘制了RDA排序图用于展示。配合plot()函数使用，除了前述scores()（提取坐标）、arrows()（绘制箭头）等命令外，还有point()（绘制散点）、text()（控制标签字体）等。通常，这几个命令相互配合使用，一个美观的RDA排序图即可呈现。
我们使用plot()，对上文中经过解释变量前向选择后所得的简约tb-RDA模型“rda_tb_forward_r”（仅使用它作为示例展示可视化方法，并没有说它是最终的模型，尽管它是通过多步计算所得，但本文中并没有结合生态学意义去评估或解释它），绘制双序图展示其中的前两个约束轴。不过在作图前，首先需要确认以下信息。
（1）根据置换检验，我们已知前两个约束轴是显著的，可以使用。
（2）根据“RsquareAdj(rda_tb_forward_r)$adj.r.squared”可知简约模型的校正后R2为“0.5814848”，然后再根据“summary(rda_tb_forward_r)”的结果可知约束轴RDA1和RDA2的承载方差分别占约束模型总方差的“0.73786”和“0.168480”，那么由此推断约束轴RDA1和RDA2的解释量分别约为“0.5814848×0.73786 = 42.91%”和“0.5814848×0.168480 = 9.80%”。解释量可观。
（3）根据所关注的重点，确定使用哪种标尺。例如这次想重点观测样方间的关系及其环境梯度，那么就以I型标尺为例展示。
（4）可首先使用“plot(rda_tb_forward_r, scaling = 1, display = c('wa', 'cn'))”查看简约作图结果，排序空间内的样方分布是否区分明显，环境变量对排序空间的贡献是否达到了预期等。
确认都没问题后，就可以作图了，一个简单的示例如下。
#plot() 作图示例，详情 ?plot.cca#以前向选择后的简约模型 rda_tb_forward_r 为例，展示前两轴，I 型标尺，双序图，默认使用物种加权和计算的样方坐标
#注意 RDA 轴的解释率，需要手动结合校正后的 R2 计算下
png('rda_test1.png', width = 600, height = 600, res = 100, units = 'px')
plot(rda_tb_forward_r, display = c('wa', 'cn'), choices = 1:2, scaling = 1, type = 'n', main = 'RDA双序图，I型标尺', xlab = 'RDA1 (42.91%)', ylab = 'RDA2 (9.80%)')
text(rda_tb_forward_r, display = 'cn', choices = 1:2, scaling = 1, col = 'blue', cex = 0.8)
points(rda_tb_forward_r, display = 'wa', choices = 1:2, scaling = 1, pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
dev.off()
其中，plot()控制全局；text()定义了环境变量向量的属性；points()定义了样方坐标散点的属性，由于示例数据中27个样方共分为3组，我们使用3种颜色表示出。
ggplot2作图
对于习惯ggplot2做图方式的同学们来说，我们只需要将样方排序坐标、物种排序坐标以及环境因子的排序坐标分别提取出构建数据集后，即可使用ggplot2绘制排序图了。
#以前向选择后的简约模型 rda_tb_forward_r 为例，展示前两轴，I 型标尺，双序图，默认使用物种加权和计算的样方坐标
#提取样方和环境因子排序坐标，前两轴，I 型标尺
rda_tb_forward_r.scaling1 <- summary(rda_tb_forward_r, scaling = 1)
rda_tb_forward_r.site <- data.frame(rda_tb_forward_r.scaling1$sites)[1:2]
rda_tb_forward_r.env <- data.frame(rda_tb_forward_r.scaling1$biplot)[1:2]
#读取样本分组数据（附件“group.txt”）
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#合并样本分组信息，构建 ggplot2 作图数据集
rda_tb_forward_r.site$sample <- rownames(rda_tb_forward_r.site)
rda_tb_forward_r.site <- merge(rda_tb_forward_r.site, group, by = 'sample')
rda_tb_forward_r.env$sample <- rownames(rda_tb_forward_r.env)
#ggplot2 作图
library(ggplot2)
p <- ggplot(rda_tb_forward_r.site, aes(RDA1, RDA2)) +
geom_point(aes(color = group)) +
scale_color_manual(values = c('red', 'orange', 'green3')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), plot.title = element_text(hjust = 0.5), legend.key = element_rect(fill = 'transparent')) +
labs(x = 'RDA1 (42.91%)', y = 'RDA2 (9.80%)', title = 'RDA双序图，I型标尺', color = '') +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_segment(data = rda_tb_forward_r.env, aes(x = 0, y = 0, xend = RDA1,yend = RDA2), arrow = arrow(length = unit(0.1, 'cm')), size = 0.3, color = 'blue') +
geom_text(data = rda_tb_forward_r.env, aes(RDA1 * 1.1, RDA2 * 1.1, label = sample), color = 'blue', size = 3)
p
#ggsave('rda_test2.pdf', p, width = 5, height = 4)
ggsave('rda_test2.png', p, width = 5, height = 4)
备注：若觉得箭头长度太长或太短，影响观测，可以同比例缩放所有箭头的长度后再绘制出来，这种做法是被允许的（上文中在介绍排序图时也提到了这种做法的可行性；而且能看到上述plot()函数的默认出图，也是这样做的，函数默认使用一个内部常数对结果重新标定，这时变量向量的箭头长度不完全等于原始值，而是等比例原始值）。
对应分析（CA）和去趋势对应分析（DCA）在群落分析中的应用